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Introduction. 


The  standard  procedure  for  computing  the  moving  average  method  of 
moments  parameter  estimates  for  an  ARMA  process  involves  the  numerical 
solution  of  a  system  of  non-linear  equations.  In  this  paper  it  will  be 
shown  how  the  so-called  inverse  autocorrelations  introduced  by  Cleveland 
[3]  can  be  utilized  to  compute  the  moving  average  parameter  estimates  as 
the  solution  to  a  system  of  linear  equations,  and  that  that  method  is 
faster  and  more  accurate  than  solving  the  non- linear  equations  alluded 
to  above.  In  addition,  the  procedure  described  below  ensures  that  an 
invertible  solution  is  obtained.  It  also  provides  a  natural  procedure 
for  obtaining  estimates  when  the  usual  estimates  do  not  exist.  These 
"estimates"  are  shown  to  provide  reasonable  initializing  values  for  cal¬ 
culating  the  maximum  likelihood  parameter  estimates. 


Text. 

The  method  of  moments  (or  preliminary  estimates,  see  Box  and  Jenkins 
[2]  A. 6. 2)  may  be  defined  as  follows: 


Definition  1.  Let  p(j)  be  the  sample  autocorrelation  function  of  a 


stochastic  process  and  let  be  the  solution  to 


p(m+l)  =  ((i^p(m) 


+  •••+  (J)^p(m-n+l) 


(1) 


p(m+l)  =  (|)j^p(m+n-13  +•••+  (|)^p(m) 


Futher  let  9,,..., 9  be  the  invertible  solution  (i.e.,  the  solution  which 
1  m 


l-ez-...-©  z  has  all  of  its  roots  outside  the  unit  circle)  to  the 
mm 


system  of  equations 


(2) 
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(where  p(j  ;Oj, . . .  ,o^,  is  the  autocorrelation  function  of  ®J^icntioiT_ 

ARMA(n,m)  process  with  autoregressive  parameters  aj,...,a^  and  moving 

Dlstrlbul j  on/ 
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average  parameters  6^ ,  •  ■ .  ,6|j|)  • 
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We  call  . the  method  of  moments  parameter 

estimators  of  an  ARMA(n,m)  process  and 

K|l.0,e2^i“ . 0 

s(ai)  =  - 

1- J  -•••-  i 

1  n  ' 

is  called  the  method  of  moments  ARMA(n,m)  spectral  density 
estimator,  where  K  is  chosen  so  that 
.5  ^ 

/  s  (dj)  do)  =  1 . 

-.5 

A  ^  ^ 

Verbally  and  0j^,...,0|jj  are  those  parameter  values  whose 

A  A  A 

theoretical  autocorrelation  function  p(j  ;(J), , . . .  ,6  ,  0„...,0  ) 

X  n  X  m 

A 

agrees  with  p(j)  for  j  =  l,...,m+n. 

A 

As  is  well-known,  for  m  >  0,  even  if  p(j)  is  a  positive  defi- 

A  A  A  A 

nite  function,  a  solution  for  4>,  .  9,,..., 9  need  not 

exist  (e,g.,  forn  =  0,  m=l,take.5  <  p  Cl)  <  1).  In  that  case 
we  say  that  the  method  of  moments  parameter  estimators  and  the 
method  of  moments  spectral  density  estimator  do  not  exist. 

In  Morton  [7  ],  a  spectral  density  estimator  called  the 
G- spectral  estimator  was  introduced  (actually  it  was  called  the 
modified  G-spectral  estimator,  it  being  a  modification  of  a 
spectral  estimator  introduced  by  Gray  [ 5  ]  and  studied  by  Gray, 
Houston  and  Morgan  ( 6  ] ) , 

The  G-spectral  estimator  may  be  defined  as  follows. 

A 

Definition  2.  Let  p(j)  be  an  estimator  for  the  autocorrelation 

function  p(j)  of  a  stochastic  process.  Further  let  f^  ■ 

and  F,  ■  E  f  ,  where  k  ■  max{l,n-m-l}  and  m  and  n  are  given, 

J  v— k  ^ 


_ 
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non-negative  integers.  We  then  define 

«  2Real(e„(F  )-F  )  +  1, 
n,ni  n  m  u 

where 


F 
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for  n  >  1  and  e-(F  )  =  F  . 

—  0  m  m 

The  following  properties  of  the  spectral  estimator  introduced 
above  are  proved  in  Morton  [  7 ] . 

A  /V  A  >N  A 

Theorem  1.  Let  G_  _(a)),  s(u)),  and  (j>, , . . .  be  as 

- —  — .  n  i  m 

defined  above.  Define  j  =  l,...,n  . 

Then 


A 

(i)  if  s(u))  exists  then 


s(u)D 

(li) 


n,m  . 


2Real 


F  -otF  ••••-otF 
m*l  °1  m  •  n  m-nfl 

1  -  <*1 . 


where  the  F^  are  as  given  in  Definition  2. 

The  first  result  above  shows  that  G  _(u))  is  an  ARMA  spectral 

n,m 

estimator  which  does  not  require  calculation  of  the  moving  average 

parameter  estimates.  The  second  result  provides  a  simple  formula 

for  calculating  G  (u) .  The  calculation  formula  we  note  does 
n,m 

require  the  calculation  of  the  autoregressive  parameter  estimates 

(from  (1),  for  instance);  however  (1)  is  a  linear  system  of 

equations  which  may  be  solved  rather  easily.  Calculating  G  (uj)  by 

rt  f  in 

(ii)  provides  a  much  more  efficient  means  of  calculating  G  (w) 

n,m 

than  does  direct  evaluation  of  the  above  determinants  at  every 
desired  frequency. 
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We  are  now  in  the  position  to  give  the  result  we  will  utilize 

in  the  calculation  of  . 

1  m 

A  A  Ak 

Theorem  2.  Let  (fc,  ,  9,  .....6  and  G  fu)!  exist  as  defined 

-  1  n  1  m  n,m 

above.  If  G  (oi)  >  0  for  all  40  and 
n,m 


ci(k)  =  / 

0 


.5 


cos(2Triu)k) 


do) 


then  0,,...,0  is  the  solution  to 
1  m 

ci(l)  =  9jCi(0)  +•••+  0^  ci(l-m) 

ci(m)  =  0,ci(m-l)  +•••♦  0  ci(0)  . 
i  m 


(3) 


Proof.  From  the  results  above 


- - ^ -  =  — - - ^ (say). 

We  note  that  f(a))  is  proportional  to  the  spectrum  of  a  AR(m)  process 
with  parameters  9j^,...,9^.  Thus, 


ci(j)  =  /  cos(2iTiu))  f(u))  du  =  y  /  e  dw  , 

0  -.5 

and  we  note  that  ci(j)  is  proportional  to  the  autocorrelation  func- 

Ak  A^ 

tion  ofanAR(ra)  process  with  parameter  values  0j,...,9^.  So  it 

A, 

follows  that  9,,...,0  is  the  solution  to  (3). 

1  m 

In  this  paper,  we  will  refer  to  algorithms  based  on  solving 
the  system  of  equations  in  (3)  as  Method  1.  We  note  from  the  above 
Theorem  that  Method  1  will  always  yield  an  invertible  moving 
average  operator  (an  algebraic  proof  that  the  operator  resulting 
from  (3)  is  invertible  is  given  in  Morton  [7]). 
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The  standard  procedure  for  calculating  the  moving  average 
parameter  estimates  is  via  an  algorithm  based  on  a  solution 
to  (2) .  The  standard  procedure  for  that  calculation  requires 
an  iterative  solution  to  a  nonlinear  system  of  equations  Csee 
Box  and  Jenkins  [2  ],  pp.  201  ff) .  Algorithms  based  on  finding 
the  solution  to  (2),  we  will  call  Method  2.  The  two  methods 
may  be  compared  by  the  computational  speed  and  accuracy  which 
they  afford  and  as  to  whether  or  not  an  invertible  solution 
is  guaranteed.  First  we  consider  the  question  of  invertibility . 

It  is  well-known  that,  even  if  a  solution  to  (2)  exists, 
it  is  not  unique.  (See  Box  and  Jenkins  [2  ],  pp.  198-199.) 
Uniqueness  is  then  guaranteed,  typically,  by  requiring  that 
the  polynomial 

A  A  A 

9(z)  *  1  -  0,2  -•••-  6  2 
1  mm 

A 

have  all  of  its  roots  outside  the  unit  circle  (i.e.,  9(B)  is 
required  to  be  invertible) .  Thus  a  numerical  algorithm  which 
merely  requires  that  the  equations  in  (2)  be  satisfied  does  not 
guarantee  an  invertible  moving  average  operator.  In  particular, 
one  commonly  used  subroutine,  FTMPS  in  IMSL,  does  not  guarantee 
an  invertible  solution.  As  noted  above.  Method  1  always  yields 
an  invertible  solution. 

To  compare  the  computational  speed  and  accuracy  of  Method 
1  and  Method  2,  we  consider  a  comparison  between  the  IMSL  sub¬ 
routine  FTMPS  and  a  program  written  by  Morton  which  calcu¬ 
lates  the  ci(j)  by  Simpson's  rule  with  Richardson  extrapolation 
(see,  for  instance,  Dahlquist  and  Bjork  [4],  pp.  269  ff) . 

To  make  the  comparison,  the  true  autocorrelations  were 
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input  for  moving  average  orders  q  =  2,3,4,5,6  and  the  two  methods 
were  compared  for  accuracy  of  the  calculated  coefficients  and 
for  the  CPU  time  required  to  make  the  calculation.  Method  1  was 
applied  using  grid  sizes  of  20,40,60,80,100  in  the  numerical  inte¬ 
gration  scheme.  A  summary  of  the  results  is  given  in  Table  1. 

Note  that,  for  these  models.  Method  1  is  faster  and  more  accurate  in 
every  case  except  for  q  =  3  with  a  grid  size  of  100  for  which  Method 


2  was  slightly  faster  (though  much  less  accurate) .  For  the  q  =  5 
case,  we  note  that  Method  2  failed  to  yield  an  invertible 
operator,  and,  for  the  q  =  4  case,  we  note  that  Method  2  failed 
to  converge  to  a  solution  at  all. 

/\  /v 

In  the  simulations  above,  a  solution  for  0,,...,0  exists. 

1  m 

However,  it  is  well-known  that  that  need  not  be  the  case.  In 
fact,  for  the  pure  moving  average  case,  a  solution  exists  if, 
and  only  if, 

m  /V 

f(a))  =  1  +  2  p(j)  cos(2ira)j)  >0,  -.5  £  uj  £  .5 

(see,  for  instance,  0.  D.  Anderson  [i  ]  pp.  137  ff . ;  for  mixed 
processes,  let  f(u))  =  G  and  add  the  condition 

T1 1  *11 

A 

that  $(B)  be  a  stationary  operator).  Thus,  if  f(u>)  <  0  for  some 
u),  no  method  of  moments  solution  exists.  One  solution  to  that 
difficulty  is  to  replace  f(a>)  by 

A 

^  -  1  +  2  E  cos(2TTaij) 

j=l 


g(u) 


1  +  c 


where  c  is  chosen  so  that  g((o)  >  0. 

In  the  m=l  case,  for  instance,  an  invertible  solution 

^  A 

exists  if,  and  only  if,  -.5  <  p(l)  <  .5.  Then,  if  p(l)  is 
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is  outside  that  interval,  the  adjustment  above  simply  "shrinks"  p(l) 

A 

to  P  ,  where  c  is  chosen  so  that  p(lV(l  +  c)  is  in  the 
1  c 

admissible  region. 

The  procedure  described  above  for  modifying  the  estimated 
autocorrelations  in  order  to  obtain  a  solution  is  not  inherently 
restricted  to  either  Method  1  or  Method  2.  However,  it  is  not 
utilized  by  any  of  the  commonly  used  algorithms  which  employ 
Method  2  (in  particular,  it  is  not  utilized  by  FTMPS  in  IMSL) ;  so 
its  use  will  only  be  considered  here  in  conjunction  with  Method  1. 

For  the  remainder  of  this  paper.  Method  1  will  refer  to  calculating 

yN  A 

®l'’‘”®m  equations  in  (3)  where  the  function 

fCoj)  = 

is  shifted  upward  if  it  takes  on  non-positive  values.  The  precise 
statistical  properties  of  the  above  estimator  are  unknown.  However, 
as  is  illustrated  below,  it  often  provides  "reasonable"  initializing 
values  for  a  maximum  likelihood  estimation  procedure. 

A  A 

For  many  practitioners,  the  calculation  of  is  per¬ 

formed  solely  for  the  purpose  of  providing  initializing  values 
for  an  iteratively  calculated,  but  more  statistically  efficient 
estimation  routine  (e.g.,  maximum  likelihood).  We,  thus, 
now  consider  the  effect  of  the  above  results  on  the  computation 
time  and  the  number  of  iterations  required  to  calculate  the 
maximum  likelihood  estimate  using  Method  1  as  against  the  standard 
procedure  of  using  Method  2  to  calculate  initializing  values. 

To  that  end,  5  realizations  of  length  200  were  simulated  from  , 
each  of  the  5  models  used  above.  Maximum  likelihood  estimates 
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were  calculated  for  each  realization  using  both  Method  1  and 
Method  2  to  obtain  the  initializing  values.  Both  the  number 
of  iterations  required  for  convergence  of  the  process  and  the 
total  CPU  time  for  the  5  realizations  from  each  separate  model 
were  recorded. 

To  summarize  the  earlier  results,  we  record  the  ways  in 
which  the  two  methods  differ: 

(1)  Method  1  is  typically  both  faster  and  more 
computationally  accurate. 

(2)  Method  1  always  yields  a  solution  and  Method 
2  need  not.  Further  it  requires  Method  2 
longer  to  determine  that  no  solution  can 

be  obtained  than  it  does  for  Method  1  to 
adjust  the  equations  and  deterndne  a 
solution  to  the  adjusted  set  of  equations. 

(3)  Method  1  ensures  an  invertible  solution, 
while  Method  2  does  not. 

We  summarize  the  results  of  the  simulation  in  Table  2,  In 
every  case,  and,  for  each  partition  size  used,  initializing  by 
Method  1  required  less  computation  time  than  initializing  by 
Method  2.  The  difference  is  primarily  due  to  the  fact  that  Method 
1  is  faster  computationally  and  always  gives  a  solution.  We  also 
note  that  in  the  case  in  which  Method  2  yielded  a  non-invertible 
operator,  the  maximum  likelihood  procedure  required  many  more 
iterations  than  the  procedure  required  when  initialized  by  an 
invertible  operator. 
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In  sunanary,  then,  we  have  introduced  an  alternative  to  the 
standard  procedure  for  calculating  the  moving  average  parameter 
estimates.  This  procedure  is  computationally  faster  and  more 
accurate  than  the  standard  procedure.  It  also  ensures  a  solution, 
as  well  as  ensuring  that  that  solution  be  invertible.  That  was 
shown  to  be  important,  in  speeding  up  the  convergence  of  the  max¬ 
imum  likelihood  estimates  calculation  when  used  as  initializing 
values  for  the  iterative  calculation  of  those  estimates. 


ES’TV-'U: 
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TABLE  1 

Approximated  Values  of  9^  Using  Various  Partition  Sizes  in  Simpson's  Rule 
Compared  to  Standard  Methods  for  Various  MA(q)  Processes 


q  =  2 


True  Parameters (9 . ) 

1 

Method  1 

Partitions 

1 . 3 

-.42 

• 

CPU  TIME 

20 

1.29977 

-.41950 

.014 

40 

1.29996 

-.41993 

.016 

60 

1.30000 

-.42000 

.019 

80 

1 . 30000 

-.42000 

.022 

100 

1.30000 

-.42000 

.026 

Method  2 

1.28450 

q  = 

-.41493 

3 

.032 

True  Parameters 

Method  1 

Partitions 

1.38 

-1.15 

.6 

20 

1.37519 

-1.14078 

.59026 

.019 

40 

1.38013 

-1.15032 

.60027 

.024 

60 

1.38006 

-1.15012 

.60013 

.025 

80 

1.38001 

-1.15002 

.60005 

.031 

100 

1.38000 

-1.15000 

.60000 

.038 

Method  2 

1.38109 

-1.14937 

.60074 

.054 

q  = 

4 

True  Parameters 

Method  1 

Partitions 

-.95 

-.9 

- .  35 

20 

-.94726 

.89509 

-.84354 

-.79266  .019 

40 

-.94726 

-.89428 

-.84105 

-.73760  .026 

60 

-.94922 

-.89825 

-.84708 

-.79578  .035 

80 

-.94986 

-.89964 

-.84935 

-.79893  .037 

100 

-.95001 

-.90001 

-.84997 

-.79991  .045 

Method  2  No  Convergence  .081 


TABLE  1  (Con't.) 
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q  a  5 


True  Parameters 

Method  1 

-1.85 

-1.75 

-1.66 

-1.56 

-.72 

CPU 

Time 

Partitions 

20 

-1.840 

-1.733 

-1.640 

-1.536 

-.703 

.022 

40 

-1.846 

-1.740 

-1.644 

-1.538 

-.707 

.030 

60 

-1.849 

-1.748 

-1.655 

-1.553 

-.716 

.041 

so 

-1.8500 

-1.7497 

-1.6591 

-1.5586 

-.7192 

.043 

100 

-1.8501 

-1.7501 

-1.6600 

-1.5600 

-.7200 

.053 

Method  2 

-2.031 

-1.917 

-1.822 

-1.711 

-.862* 

.071 

True  Parameters 

Method  1 

-1 

-.18 

-.17 

-.15 

.61 

.61 

Partitions 

20 

-.995 

-.183 

-.  175 

-.155 

.594 

.590 

.023 

40 

-.997 

-.177 

-.166 

-.145 

.600 

.597 

.034 

60 

-.996 

-.1792 

-.1684 

-.1479 

.6072 

.6062 

.044 

80 

-1.0000 

-.1799 

-.1696 

-.1494 

.6094 

.6092 

.050 

100 

-1.0001 

-.1300 

-.1699 

-.1499 

.6100 

.6100 

.059 

Method  2 

-.9995 

-.1791 

-.1712 

-.1489 

.6094 

.6101 

.082 

* 

This  operator  is  non-invertible. 


TABLE  2 


12 


Table  of  the 

q  »  2 

Realization 

Total  CPU 
Time 

q  =  3 

Realization 

Total  CPU 
Time 

q  a  4 

Realization 

Total  CPU 


Number  of  Iterations  Required  for  Method  1  and  Method  2 


Method  1 


Method  2 


^  PartTtions 

20 

40 

100 

200 

9 

9 

9 

9 

9 

11 

11 

11 

11 

24* 

13 

13 

13 

13 

17* 

9 

9 

9 

9 

9 

13 

13 

13  • 

13 

17* 

2.775 

2.780 

2.859 

2.935 

4.149 

Method 

1  _ 

Method  2 

Partitions 

20 

40 

100 

200 

17 

17 

17 

17 

17* 

21 

25 

22 

21 

19* 

11 

11 

11  . 

11 

11 

15 

13 

13 

13 

15* 

15 

15 

15 

15 

15 

5.552 

5.659 

5.558 

5.698 

6.177 

Method  1 


Method  2 


20  40 _ 1^0 _ 200 


1 

23 

21 

19 

19 

53** 

2 

15 

15 

15 

15 

15 

3 

9 

9 

9 

9 

13* 

4 

15 

19 

17 

17 

19* 

5 

23 

23 

23 

23 

25* 

7.940 

8.279 

7.961 

8.143 

12.714 

Time 
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TABLE  2  (Cont'd.) 


q  .  5 


Realization 


Total  CPU 
Time 


Method  1 

Method  2 

20 

Partitions 

40  100 

200 

1 

63 

53 

35 

33 

73* 

2 

45 

45 

57 

43 

41* 

3 

33 

35 

35 

35 

45* 

4 

23 

35 

21 

21 

43* 

5 

41 

35 

3S 

35 

37* 

24.296 

24.290 

22.163 

20.454 

31.477 

q  =  6 


Method  1 
Partitions 


Method  2 


20 

40 

100 

200 

1 

13 

17 

17 

17 

15* 

2 

13 

13 

11 

13 

11* 

Realization 

3 

17 

13 

13 

13 

13* 

4 

25 

21 

21 

21 

15* 

5 

13 

7 

11 

11 

15* 

Total  CPU 

12.905 

11.610 

11.521 

12.500 

16.689 

Time 


Summar\'  of  the  efficiencies  for  choosing  starting  values  for  the  maximum 
likelihood  estimation  routine  of  Method  1  and  Method  2.  The  values  given 
in  the  Table  are  the  number  of  iterations  required  for  convergence  to  a 
solution. 

* 

No  solution  exists  to  (2.2). 

Hr* 

A  non-invertible  solution  to  (2.2)  was  obtained. 
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